First-principles study of the lattice and electronic structures of TbMn205 
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The structural, electronic and lattice dielectric properties of multiferroic TbMn205 are investi- 
gated using density functional theory within the generalized gradient approximation (GGA). We use 
coUinear spin approximations and ignore the spin-orbit coupling. The calculated structural param- 
eters are in excellent agreement with the experiments. We confirm that the ground state structure 
of TbMn205 is of space group Pb2im, allowing polarizations along the 6-axis. The spontaneous 
electric polarization is calculated to be 1187 nCcm~^. The calculated zone-center optical phonons 
frequencies and the oscillator strengths of IR phonons agree very well with the experimental val- 
ues. We then derive an effective Hamiltonian to explain the magnetically-induced ferroelectricity in 
this compound. Our results strongly suggest that the ferroelectricity in TbMn2 05 is driven by the 
magnetic ordering that breaks the the inversion symmetry, without invoking the spin-orbit coupling. 

PACS numbers: 75.25.-l-z, 77.80.-e, 63.20.-e 
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I. INTRODUCTION 

Recently, a large class of manganese oxides (RMnOa jii^ 
and RMnaOs^iii^ with R=Y, Tb, Dy, etc.) has been dis- 
covered to be multiferroic, with strong magnetoelectric 
(ME) coupling. The ME coupling leads to various novel 
physical effects, such as the "colossal magnetodielec- 
tric" (CMD) effects and magneto-spin-flop effects.— i^i^ 
For example, in TbMn2 05,'^'^ the ME coupling is so 
strong, that the electric polarization can be reversed 
by applying a magnetic fields The remarkable ME ef- 
fects revealed in these materials have attracted great 
attention^'^'^'^'^'^'^ d°did^ because of the fascinating 
physics and their potential applications in novel multi- 
functional ME devices. 

The crystal structure of TbMn2 05 is orthorhombic, 
with four chemical formula units per primitive cell (32 
atoms in total), containing Mn^'+Og octahedra and 
Mn'^+Os pyramids, as shown in Fig[TJ TbMn205 shows 
several magnetic phase transitions accompanied with the 
appearance of electric polarizations and dielectric anoma- 
lies, when cooling down from the room temperaturei^'^ 
Starting from an incommensurate antimagnetic (ICM) 
ordering at Tn — 43K with a propagation vector k 
(~0.50, 0, 0.30), the structure locks into commensurate 
antimagnetic (CM) state at Tcm — 33K with k=(0.5, 
0, 0.25), during which spontaneous polarization arises at 
TpE — 38K4ii. When the temperature lowers to Tjcm = 
24 K, magnetic order becomes ICM again, with a sud- 
den decrease of polarization and a jump of the k vector 
to (0.48, 0, 0.32). The spontaneous polarization increase 
again, when continuing to cool to about lOKi^ During 
the magnetic phase transitions, a peak at TpE and a 
step at TicM of the dielectric constant was observed;^ 
indicating strong ME coupling in this compound. It was 
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demonstrated the electric polarization can be reversed by 
applying magnetic fieldi^ 

Experimental data show that the structure of 
TbMn205 has space group Pbam^ which includes spa- 
tial inversion (i?~^) symmetry. It is therefore puzzling 
that the material can develop spontaneous electric polar- 
izations. It has been suspected^i^ that the actual symme- 
try group of TbMn205 is Pb2im, allowing polarization 
along the b axis. Indeed, there are several experiments 
supporting this hypnosis. '*'-'-'^ For example, some Raman 
modes were found to be IR active in TbMn2 05,^° and 
the anomalies of atomic displacement parameters (ADP) 
have been observed.* Nevertheless, no direct evidence of 
the lower symmetry has yet been found.— "^ 

Theoretically, the microscopic origin of the strong ME 
coupling and the electric polarization is still under in- 
tensive debates.-^i^ii^ The ME coupling could originate 
either from the symmetric superexchange interactions, 
or from the antisymmetric exchange interactions.- The 
antisymmetric exchange comes from the spin-orbit cou- 
pling, and the noncoUinearity of the spin structure is an 
essential ingredient for this mechanism4 d^d^i^^ However 
it was shown in Ref. 4 that the largest electric polariza- 
tion in TbMn205 is associated with the CM state that is 
almost coUinear J^ In our recent work)i& we determined 
the ground-state structure of TbMn205 using the first- 
principles methods. The results show that the ground- 
state structure is indeed of polar space group Pb2im and 
the electric polarization equals 1187 nC-cm^^. In the cal- 
culations, we use coUinear spin approximation and ignore 
the spin-orbit interaction, suggesting that ME coupling 
in TbMn205 is due to the symmetric superexchange in- 
teractions. 

The aim of the present work to examine rigorously the 
ground state structural, electronic and lattice dielectric 
properties of TbMn205 using first-principles calculations 
to provide solid ground for further investigations. The 
rest of paper is organized as follows. After a brief dis- 
cussion of the first-principles methods and the approxi- 



TABLE I: The calculated lattice constants compared with 
experimental data. FM and L are the structures of spin 
configurations a and g in Fig. (2] respectively. 



lattice const. (A). 



FM 



L 



Exp. 



a 
b 
c 



7.3170 


7.3014 


7.3251 


8.5269 


8.5393 


8.5168 


5.6611 


5.6056 


5.6750 



mations used in the calculations in Sec. [TTl we provide a 
detailed analysis of the ground-state structural and elec- 
tronic properties in Sec. IIIII IIVI In Sec. |Vl we calculate 
all zone center optical phonon frequencies and the oscil- 
lator strengths of IR modes. The results are in excellent 
agreement with the known experimental IR and Raman 
spectra. In Sec. |Vll we calculate electric polarization in 
TbMn205. We then derive an effective Hamiltonian to 
explain the microscopic mechanisms of the ferroelectric- 
ity and the giant magnetoelectric coupling. We conclude 
in Sec. IVlTl 




FIG. 1: (Color online) Structure of TbMn205 unit cell, show- 
ing Mn*"'"06 Octahedra and Mn^'^Os pyramids. 



III. GROUND-STATE STRUCTURE 



II. METHODOLOGY 

Our calculations are based on the standard density- 
functional (DFT) theory with spin-polarized general- 
ized gradient approximation (GGA). We adopt Perdew- 
Burke-Ernzerhof functional"'^^ implemented in the Vi- 
enna ah initio Simulations Package (VASP).^*'^^ A 
plane-wave basis and projector augmented-wave (PAW) 
pseudopotentials^ are used, with Mn ipidAs, and Tb 
hpbdQs electrons treated self-consistently. A 500 eV 
plane-wave cutoff results in good convergence of the to- 
tal energies. We relax the structure until the changes 
of total energy in the self-consistent calculations are less 
than 10~^ eV, and the remaining forces are less than 1 
meV/A. Experimentally, TbMn2 05 is found to be incom- 
mensurate anti- ferromagnetic (AFM) below 24 K, with 
the propagation vector k k, (0.48, 0, 0.32). To accommo- 
date the magnetic structure, one needs a huge supercell, 
which is computationally prohibitive. Instead, we use a 
2x1x1 supercell, equivalent to approximating the propa- 
gation vector k — (0.5, 0, 0). The validity of this approxi- 
mation has been justified in our previous workfi^ For the 
supercell we used, a 1 x 2 x 4 Monkhorst-Pack k-points 
mesh converges very well the results. 

It was demonstrated in Ref. [j that the largest elec- 
tric polarization is associated with a commensurate mag- 
netic (CM) state that is almost collinear.^^ Therefore, 
in the calculations, we use the coUinear spin approxi- 
mation and ignore the spin-orbit coupling. Our results 
agree very well with the known experiments, indicating 
that these approximations capture the essential physics 
in TbMnsOs. 



We start the structural relaxation from the ex- 
perimental structural parameters-i^ To determine the 
ground state structure, we tried various spin config- 
urations (SCs). Eight selected SCs, four ferromag- 
netic/ferrimagnetic (FM) and four antiferromagnetic 
(AFM), are shown in Figl^l Here, we consider only the 
spins of Mn'^+ and Mn^+ ions. The total energy of each 
SCs is calculated with full relaxations of the (electronic 
and lattice) structures. The stablest SC, i.e., the SC of 
lowest total energy, labeled as g in Fig[2] is AFM and 
is identical to the spin structure proposed in RefiJ. SC 
g has an energetically degenerate SC,— labeled as g' in 
Fig. [3 In both SCs g and g' , Mtl^+ form an AFM 
square lattice in the ah plane, whereas Mn'^^ couples to 
Mn^+ either antiferromagnetically via J4 along a axis or 
with alternating sign via J3 along the h axis. Mn'^^ ions 
in two connected pyramids also couple antiferromagneti- 
cally through J5. Here, we adopt the notations J3, J4 and 
J5 from Ref. i, and define the J3 to be the Mn'*+- Mn3+ 
superexchange interaction through pyramidal base cor- 
ners, and J4 the superexchange interaction through the 
pyramidal apex, as indicated in Fig. [3] We label the two 
different Mn^+ chains along the a axis I, II, respectively, 
also following Ref. [J. The magnetic structure of g' can 
be obtained from g by shifting chain II to the right (or to 
the left) by one unit cell along the a axis.*^ The exchange 
integrals J3, J4, J5, were fitted via a Heisenberg model 
assuming nearest neighbor (NN) coupling. We find that 
all the exchange energies are of AFM type, i.e., negative 
and |J4|,|J5| ^ \J^\. Therefore, the spins must couple 
via J4 and J5 anti-ferromagnetically in the stable mag- 
netic structures. More details about the spin structure is 
given in Appendix \K\ 




FIG. 2: (Color online) Spin configurations of Mn^+, Mn*+ 
ions in the ab plane. 



The lattice structure relaxed from SC g lead to struc- 
ture L in Ref. \IQ, whereas SC g' gives the structure 
R. The calculated lattice constant^i^ of the ground state 
structure L (and R) are listed in Table [H and are in 
very good agreement with the experiments.— The er- 
rors of the lattice constants are about 1%, typical er- 



rors for GGA. The lattice constants calculated for the 
full FM configuration (SC a), are also hsted in Table HI 
which are somewhat different from those of structure L. 
The calculated Wyckoff positions (WPs) for the struc- 
ture L are given in Table [TTl comparing with the ex- 
perimental data^ whereas the structure i? is a mirror 
image of L about the ac-plane. The calculated WPs 
positions are extremely close to what was obtained ex- 
perimentally. However, the WPs of structure L split by 
small amount from the WPs of the (Pbam) experimen- 
tal structure, and lower the structural symmetry to the 
long searched Pb2im polar group. ^ The relationships be- 
tween the WPs of TbMn205 under Pbam symmetry and 
under the Pb2im symmetry are given in Table ITTTl For 
example, under the Pbam symmetry, Mn'^+ has one WP 
(notion h),with four equivalent sites. It splits into two 
WPs b(l) and b(2), each having two equivalent sites, 
when the symmetry lowers to Pb2im. 

We artifically construct a high symmetry structure H 
by symmetrizing structure L and R according to the 
Pbam symmetry. ^^ The WPs of H together with the 
atomic displacements from H to L are also given in Table 
im where 5a,, 5h and 6c denote the atomic displacements 
from the high symmetry positions along the a, b and c 
axes respectively. The displacements along the a and c 
axes are of mirror symmetry, whereas the displacements 
along the b axis are not. As seen in Table [III all cations 
move up and anions (except Oi) move down from the 
positions of high symmetry structure H, resulting in po- 
larization. However, all the atomic displacements are 
extremely small, usually are of the order of ~ 10"'* of 
the lattice constants. Therefore the low symmetry struc- 
ture can not be directly determined experimentally, and 
only the anomalies of the atomic displacement parame- 
ters (ADPs) were observed.^ The largest atomic displace- 
ments come from Mn'^+, and O3, 6y ^ 10^"^ of the lattice 
constants, in agreement with the ADP results of Reflj. 
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IV. ELECTRONIC STRUCTURE 



Figure [4] depicts the calculated the total densities of 
states (DOS) of structure L . The DOS of spin-up and 
spin-down electrons are identical as expected for an AFM 
state. From the results, the band gap is estimated to be ^ 
0.4 eV, confirming the experimental fact that TbMn2 05 
is an insulator. However, it is well known that GGA 
greatly underestimates the band gap, especially for the 3d 
compounds, therefore the real band gap might be much 
larger. We further calculated the site-projected partial 
DOS (PDOS) for Tb, Mn d electrons and O p electrons 
in the energy range of -2.5 to 2.5 eV around the Fermi 
level, shown in FigO In structure L, the Tb, Mn'^+, 
O2, O3 and O4 each has two non- equivalent sites, and 
the PDOS are shown separately. The PDOS for other 



electrons are quite small in the selected energy range and 
therefore are not shown. We see the DOS near the Fermi 
level is mainly from Mn 3d and O 2p orbitals . 

The local magnetic moments are estimated for Mn'^+ 
to be ^ 2.37 fiB, and for Mn^+ to be ~ 1.64 /i^, in 
a good agreement with the refined magnetic momentsi^ 
The ionic radii of 0.60 A and 0.66 A are used for Mn''^+ 
and Mn'^+ respectively taken from the periodic table. 
The local magnetic moments of Mn'^+ , and Mn''+ at dif- 
ferent WPs are slightly different. O ions also have small 
induced magnetic moments, due to their hybridization 
with Mn ions. The calculated local moments of Tb are 
extremely small, whereas 1.6 /is is assumed in order to 
obtain a perfect fit of the neutron diffraction data in Ref. 
4. The difference might come from that the / electrons 
are not included as valence electrons in the Tb pseudopo- 
tential. 




FIG. 3: (Color online) Spin configurations with different combinations of chains I and II. g and g' are the ground state spin 
configurations for structure L and R respectively, whereas dw and dw' are the domain walls between structure L and R. 

TABLE II: Comparison of the calculated and measured atom positions of TbMn205 . L and H are the ground state structure and 
the high-symmetry structure respectively. \Sa,\, \Sh\, and \5c\ denote the atomic displacements from H to L. The experimental 
values are taken from Ref. 1 131. 







L {Pb2im) 






H (Pbam) 




H 


-> L(10- 


-4) 


Exp. (Pbam) 


atom 


a 


b 


c 


a 


b 


c 


\Sa\ 


\5h\ 


\5c\ 


a 


b 


c 


Tbf+ 


0.1410 


0.1733 





0.1407 


0.1732 





3.0 


1.5 





0.1399 


0.1726 





Tb^+ 


0.6404 


0.3270 











3.0 


1.5 











Mn*+ 


0.0001 


0.5003 


0.2558 





0.5 


0.2558 


0.8 


2.9 








0.5 


0.2618 


Mn?+ 


0.4012 


0.3558 


0.5 


0.4014 


0.3551 


0.5 


2.2 


6.6 





0.4120 


0.3510 


0.5 


Mn^+ 


0.9016 


0.1456 


0.5 








2.2 


6.6 











Ol 


0.0008 


0.0002 


0.2709 








0.2709 


8.2 


2.3 











0.2710 


02i 


0.1646 


0.4480 





0.1647 


0.4481 





1.2 


1.2 





0.1617 


0.4463 





022 


0.6648 


0.0517 











1.2 


1.2 











03i 


0.1560 


0.4329 


0.5 


0.1565 


0.4337 


0.5 


5.3 


7.8 





0.1528 


0.4324 


0.5 


032 


0.6571 


0.0655 


0.5 








5.3 


7.8 











04i 


0.3977 


0.2077 


0.2438 


0.3968 


0.2079 


0.2430 


8.8 


2.2 


8.5 


0.3973 


0.2062 


0.2483 


042 


0.8959 


0.2919 


0.7579 








8.8 


2.2 


8.5 










TABLE III: The Wyckoff positions (WPs) for each inequiva- 
lent atom in TbMn2 05 under the space group Pbam and and 
its subgroup Pb2\m. 



Energy (eV) 



Atoms 


WP(P6am) 


WP(P62im) 


Tb^+ 


4g 


2a, 2a 


Mn*+ 


4f 


4c 


Mn3+ 


4h 


2b, 2b 


01 


4e 


4c 


02 


4g 


2a, 2a 


03 


4h 


2b, 2b 


04 


8i 


4c, 4c 



FIG. 4: The total density of sates (DOS) of TbMn205 for 
spin-up and spin-down respectively, calculated under struc- 
ture L. The dashed line represents the Fermi level. 



We also calculated the DOS and PDOS of structure H 
and find no obvious differences from structure L, suggest- 
ing that the effects of the lattice structure distortion and 
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FIG. 5: Site-projected partial density of states (PDOS) of 
Mn and Tb d electrons, and O p electrons for TbMn2 05. 
The upper panel and lower panel show the PDOS for spin up 
and spin down respectively. The PDOS of other states are 
relatively small in the selected energy range, and therefore 
not shown. The dashed lines represent the Fermi level. 



TABLE IV: The relationships of irreducible representations 
between space group Pbam and its subgroup Pb2\m. 





Ph2im 




Pbam 


irreps 


Character 


irreps 


Character 


Ai 


IR and Raman 


As 


Raman 






B2u 


IR 


Bi 


IR and Raman 


Bag 


Raman 






Bill 


IR 


B2 


IR and Raman 


Big 


Raman 






B3„ 


IR 


A2 


Raman 


B2g 


Raman 






A„ 


Silent 



not observed directly in the experiments. However, a 
strong evidence that suggests the symmetry lowering is 
that one of the Raman-active phonons become also IR- 
active at low temperature,^^ which is forbidden by the 
higher Pbam symmetry. To elucidate the origin of the 
symmetry lowering, we calculate the frequencies of all 
zone center phonons, using frozen phonon method, to 
exam if any soft modes exist in TbMn205 and the sym- 
metry properties of all phonons. Furthermore, the va- 
lidity of our theories can be examined by comparing the 
calculated phonons with experiments. 



A. symmetry analysis 

We first carry out symmetry analyses and decompose 
the symmetry- adapted modes into irreducible represen- 
tations. The Hamiltonian remains invariant under a lat- 
tice translation of the 32-atom primitive cell combined 
with spin reversal. Therefore the force-constant matrix 
has the translational symmetry of the 32-atom primitive 
celli^ We then perform the symmetry analyses on the 
32-atom primitive cell (instead of the 64-atom AFM cell 
used in the calculations), using the SMODES symmetry- 
analysis software package i^ 

For the high symmetry structure {Pbam), all the 96 
modes, including 3 acoustic modes, are decomposed into 
8 irreducible representations (irreps): 

r{Pbam) = 9Bi„ © 15B2u © l^B^^ 13Ag 

mSBig © llB2g © UBsg © 9Au . (1) 

Among them B^u, B2u and Bi„ modes are IR active, 
polarized along the a, b and c axis, respectively, whereas 
Big, i?2g, B^g and Ag modes are Raman active and v4„ 
modes are silent. All these irreps are one dimensional, 
i.e., no degenerate phonons. The three acoustic modes 
are each in one of the three IR-active representations. As 
we see, the Raman- and IR-active modes do not couple 
under the Pbam symmetry. 

For the crystal structure processes a Pb2im, symmetry, 
the 96 phonons can be decomposed into 4 irreps: 

r(P62iTO) ^ 28Ai © 20Bi © 28^2 ® 20^2 . (2) 



the symmetry lowering in structure L to the electronic 
structure are small. 



V. LATTICE DYNAMIC PROPERTIES 

The symmetry lowering and ferroelectricity in the tra- 
ditional ferroelectrics, such as BaTiOa, PbTiOa, etc. are 
driven by the condensation of zone center soft phonons. ^^ 
In TbMn205, the low symmetry {Pb2im) structure was 



We found all modes are Raman active, among them Ai, 
Bi and B2 are also IR active, which are polarized along 
the b, c, and a axis respectively^^ The relationships be- 
tween the irreps of Pbam symmetry and its subgroup 
Pb2im, symmetry are given in Table HVl together with 
the symmetry character for each irrep. The Ai modes 
are coupled from the Ag and B211 modes and the A2 
modes are composed of the i?2g and Au modes. The 
Bi modes are the coupled B^g and i?i„ modes whereas 
the B2 modes are the coupled Big and B^u modes. 



TABLE V; Calculated phonon frequencies u), effective charges Z^ and oscillator strength 5a of A\ modes, classified according 
to their major symmetry characters and compared with the experimental data. The IR-active phonons are extracted from Ref. 
[lOl at 7K and Raman-active phonons are extracted from Ref. Y2A . 



GGA Exp. 



B2u(IR, h- 
GGA Exp 



polarized) 

Sa 
. GGA 



A, 



Raman% 



Exp. 



GGA Exp. 



(Raman] 

Zl 
GGA 



GGA 



IR% 



100.7 


97.2 


0.09 


0.17 


0.11 


0.42 


0.08 


110.1 




0.020 


0.005 


0.02 


158.0 


168.9 


0.28 


0.31 


0.44 


0.46 


0.12 


136.9 




-0.003 


~0 


0.17 


162.8 


171.9 


0.31 


0.25 


0.49 


0.30 


0.04 


221.6 


215 


0.058 


0.009 


0.94 


224.8 


222.2 


0.30 


0.20 


0.24 


0.11 


1.85 


235.1 


221 


0.066 


0.011 


1.38 


267.3 


256.8 


0.25 


0.28 


0.12 


0.17 


0.36 


312.5 


334 


0.125 


0.022 


0.97 


316.7 


333.4 


0.78 


0.37 


0.82 


0.17 


1.02 


340.2 


350 


0.094 


0.011 


2.96 


351.3 


386 


0.28 


0.15 


0.09 


0.02 


2.94 


405.6 


412 


-0.028 


~ 


0.24 


412.5 


422.3 


0.39 


0.60 


0.13 


0.28 


0.23 


445.1 




0.423 


0.126 


2.24 


439.5 


453.2 


2.58 


2.25 


4.81 


3.43 


2.65 


489.2 


500 


0.151 


0.013 


6.51 


471.0 


481.8 


1.40 


2.18 


1.23 


2.86 


5.89 


529.2 


537 


-0.104 


0.005 


6.16 


533.5 


538.2 


0.52 


0.72 


0.27 


0.25 


5.95 


612.0 


621 


-0.002 


~0 


3.28 


549.3 


567.3 


0.50 


1.09 


0.12 


0.52 


1.33 


613.5 


631 


0.030 


~0 


9.59 


625.0 


636.6 


1.01 


0.88 


0.36 


0.27 


11.20 


673.6 


693 " 


-0.107 


0.004 


8.04 


667.2 


688.3 


0.17 


0.10 


0.01 


0.003 


8.83 













"IR-active and 703 cm ^ measured by Aguilaa 



TABLE VL Calculated phonon frequencies lo, effective charges Za and oscillator strength Sx of B2 modes, classified according 
to their major symmetry characters and compared with the experimental data. The IR-active phonons are extracted from Ref. 
Hoi at 7K and Raman-active phonons extracted from Ref. |24| . 



GGA Exp. 



B3u(IR, a- 
Zl 
GGA Exp 



polarized) 
Sa 
. GGA Exp. 



B 



Raman% 



19 



GGA Exp. 



(Raman) 
GGA GGA 



IR% 



111.8 


111.9 


0.02 


0.23 


0.004 


0.59 


0.02 


148.2 


148 


0.007 


~0 


0.54 


152.2 


157.5 


0.17 


0.37 


0.18 


0.81 


0.50 


172.9 


172 


0.006 


~0 


0.13 


159.8 


164.2 


0.76 


0.57 


3.11 


1.68 


0.12 


222.1 


208 


0.054 


0.008 


10.5 


217.4 


218.5 


0.05 


0.32 


0.01 


0.30 


10.99 


243.0 


237 


0.067 


0.011 


0.97 


264.5 


254.8 


0.30 


0.94 


0.18 


1.88 


0.62 


300.7 


326 


0.007 


~0 


0.18 


317.0 


333.1 


0.03 


0.27 


0.002 


0.09 


0.15 


373.0 




0.048 


0.002 


9.50 


353.1 


364.9 


1.04 


1.39 


1.21 


2.02 


1.25 


410.4 


416 


0.019 


~0 


0.31 


379.1 


397.6 


0.56 


0.66 


0.30 


0.38 


8.44 


470.0 


485 


0.360 


0.081 


36.8 


465.2 




0.21 




0.03 




1.20 


518.3 




0.115 


0.007 


2.82 


474.4 


494.8 


0.24 


0.89 


0.03 


0.45 


37.76 


545.8 


538 


0.054 


0.001 


0.58 


509.2 




1.18 




0.75 




0.59 


587.1 




0.178 


0.013 


1.91 


595.4 


613.5 


0.90 


1.38 


0.32 


0.71 


1.28 


643.8 


673 


0.445 


0.067 


5.35 


617.1 


627.5 


1.23 


0.81 


0.55 


0.23 


4.57 


678.7 


697 


0.027 


~0 


10.0 


672.7 


704.2 


0.33 


0.42 


0.03 


0.05 


12.12 













B. Zone Center Optical Phonons 

The zone-center phonons frequencies are calculated via 
a frozen phonon technique. ^^ For each irrep T in Eq. ([1]), 
(p|) , we construct an N-p x Ny dynamical matrix from a se- 
ries of frozen-phonon calculations in which the structure 
is distorted according to each of the iVr symmetry modes 
in this irrep. For example, A^r=28 for the Ai modes in 
Eq. ([2]). All calculations were performed in the 64- atom 
AFM cell respect to the correct AFM reference ground 
state. After calculating the residual Hellmann-Feynman 
forces we displace the ions according to the symmetry 



coordinate for each mode by 0.1% of lattice constant and 
calculate the forces F" on ion i in Cartesian direction 
a due to displacements u^ of ion j in direction /3. The 
force constant matrix 



$ 



Q/3 



dFT 



(3) 



is then obtained by finite differences to the residual forces 
in the equilibrium structure to eliminate the numerical 
errors due to the residual forces. The normal modes U;^ 
and their frequencies uj\ are then obtained through solu- 



tion of the eigenvalue equation 

^ ■ ux = LufM ■ ux. (4) 

Here M"^^ = (mi/mo) 5ij Sap is the dimensionless diag- 
onal mass matrix, where rrii is the mass of atom i and 
Too is a reference mass chosen here to be 1 amu, and the 
eigenvectors are normalized according to u^-M-u^ = (5^^. 

The phonon frequencies calculated from high structure 
H are extremely close to those calculated from low sym- 
metry structure L. No soft mode is found in both high 
and low symmetry structures, excluding the possibility 
of ferroelectricity driven by the soft modes. We therefore 
show only the phonon frequencies for the ground state 
structure L only. The calculated Ai, B2, Bi phonon fre- 
quencies are given in Table El Table IVI[ and Table IVIII 
respectively, whereas the A2 phonons are given in Ta- 
ble Ivini To make a good contact with experiments, we 
divide the phonons of each irrep into two presentations 
according to their major symmetry characters (see Table 
HV)) . We found that for most of the modes the coupling 
between the major and minor characters (IR and Ra- 
man or Raman and silent) are extremely small, and such 
classification is unambiguous. All calculated phonon fre- 
quencies are in excellent agreement with the experiments 
especially for the low-frequency phonons. The errors are 
less than 5%. 

In order to check the influences on phonons from dif- 
ferent magnetic orderings and magnetic-ordering-induced 
lattice modulations, we also calculated the phonon fre- 
quencies for the FM structure, we find the differences 
are also small. 

In the following sections, we compare the calculated 
IR and Raman spectra with the experiments in details. 



1. IR spectra and dielectric properties 

Aguilar et al. measured the frequencies of IR-active 
phonon polarized along the b and a axes.^^ These phonon 
modes are compared with the calculated modes in Table 
IVl and Table |Vll respectively. 

For the IR phonons, we can further calculate their con- 
tribution to the dielectric constant. The dielectric func- 
tion tensor in the frequency range near the phonon fre- 
quencies could be written as sum 



ea/3(w) — (eoo)a/3 + {<^ph)afi{^) 



(5) 



of electronic contribution and lattice contribution, where 
a, (3 = a, b, c are the axis indices. For TbMn205, 
(eoo)aa and (eoo)b6 were measured to be 5.31 and 6.82 
respectively.— In TbMn205, the phonons are polarized 
along the a,b,c axis, and eapiuj) is diagonal. We could 
therefore ignore the subscripts a, (3, but only give the 
polarization direction. The dielectric contribution from 
the a-polarized IR-active phonons is then 



where u\ and Z^ are the mode frequencies and mode 
dynamical charges respectively, and flo^ = 'iire^/moV is 
a characteristic frequency having the interpretation of a 
plasma frequency of a gas of objects of mass to,o=1 amu, 
charge e, and density V~^ {V is the 32-atom primitive 
cell volume) i2^ The sum is restricted to the a-polarized 
phonons. The mode dynamical charges are defined as 



z: 



V 



dPg 

dux 



(7) 



eS^M^f^^E- 



zf 



^i- 



UJ^ 



(6) 



in which P^ is the polarization in the a direction due to a 
small frozen-ion amplitude of the A-th phonon modesi^ 
In practice we compute the mode dynamical charges by 
finite differences method following Ref. Q^ Once we 
have the mode effective charges, we calculate the oscilla- 
tor strengths of each IR-active modes Sx — ^o^x^/^x^ 
and the total lattice contribution to the static {uj -^ 0) 
dielectric constant eph(O) = ^oJ2x ^x^/^x = X^a '^^• 

The results of the Born effective charges and oscillator 
strengths for the IR-active modes are also listed in Table 
EllVllandlVlTl for the Ai (6-polarized), B2(a-polarized) 
and Bi (c-polarized) phonons respectively. 

Let us first look at the Ai (6-polarized) modes pre- 
sented in Table FVl There are 28 Ai modes in total, and 
14 (plus one acoustic mode) of them that are mainly i?2u 
IR modes. The calculated phonon frequencies and oscil- 
lator strengths are in very good agreement with experi- 
ments. The calculated total oscillator strength is about 
9.24 compared with the measures value 9.12. The rest 
13 modes are mainly Ag Raman modes. Under the low 
Pb2im symmetry, these modes also acquire some small 
oscillator strengths. More importantly, the 693 cm~^ Ag 
Raman mode (703 cm~^ measured from IR spectraJi), 
that was found also IR active with oscillator strength 
S'a=0.001 in the experiment, ^"^ is well reproduced in the 
calculations, with oscillator strength S'a=0.004, there- 
fore, confirming that the ground state structure is indeed 
of Pb2im symmetry.-iS, The oscillator strengths of other 
Ag modes might be covered by the adjacent i?2« modes, 
and therefore are not observed experimentally. To see 
more clearly between the coupling of the IR and Raman 
modes, we also list the percentage of the characters in 
the Table. We see the coupling between IR and Raman 
modes is fairly small. 

The results of 28 B2 phonons are given in Table IVIl 
classified into two presentation according to their major 
symmetry characters, including 14 IR (B^u) modes, 13 
(Big) Raman modes and one acoustic mode. Two of the 
-63,1 modes were missing in the experimentSf^ possibly 
due to their small oscillator strengths. Here we tenta- 
tively assign them as the 465.2 and 509.2 cm^^ modes. 
Again the phonon frequencies agree very well with the 
experiments. However, the IR oscillator strengths are 
in less agreement with the experimental values. Some 
modes have large discrepancy in the oscillator strengths 
compared to the experiments, e.g, the 264 cm^^ mode. 
Nevertheless, the calculated total oscillator strength for 
a-polarized phonons is about 6.87, in a reasonable agree- 



TABLE VII: Calculated phonon frequencies uj, effective charges Z^ and oscillator strength S\ of Bi modes, classified according 
to their major symmetry characters. The Raman-active phonons are extracted from Ref. |2J. 



Biu(IR, 


c-polarized) 




B; 


ig (Raman) 




U) 


zt 


Sa 


Raman% 


LO 




zi 


Sx 


IR% 


GGA 


GGA 


GGA 




GGA 


Exp. 


GGA 


GGA 




142.2 


0.49 


1.65 


~0 


102.6 




0.005 


~0 


0.004 


249.5 


1.01 


2.27 


0.62 


192.1 


195 


0.007 


~0 


0.03 


300.2 


0.18 


0.05 


3.28 


241.6 


244 


0.100 


0.024 


0.41 


407.4 


1.47 


1.828 


24.58 


265.5 




0.040 


0.003 


0.69 


433.9 


0.84 


0.53 


27.42 


289.0 


299 


0.107 


0.019 


1.50 


495.8 


1.20 


0.82 


36.47 


325.2 




0.098 


0.013 


1.27 


527.2 


0.39 


0.08 


2.80 


426.4 




2.25 


3.89 


48.8 


605.1 


0.74 


0.21 


0.56 


447.2 
470.0 


442 


0.116 
0.124 


0.009 
0.010 


0.45 
5.92 










505.8 


540 


0.909 


0.450 


34.5 










553.9 


577 


0.012 


~0 


2.07 



ment with the experiments value about 9.191^^ Similar 
to the Ag modes, the Big Raman modes also have small 
IR oscillator strengths due to the Ph2im symmetry. For 
most modes, the coupling between the IR and Raman 
characters are small, except two modes: the 474 cm~^ 
mode in the i?3„ presentation and the 470 mode in the 
Big presentation. This is probably because they are ac- 
cidentally degenerate in phonon frequencies, leading to 
large character mixing. 

The phonon frequencies and oscillator strengths of the 
c-axis polarized Bi phonons are given in Table IVIH in- 
cluding 8 (plus one acoustic mode) Biu IR modes and 11 
B^g Raman modes. Interestingly, there are some modes 
that have very large IR and Raman modes mixing in 
this representation. For example, the 426 cm~^ Raman 
modes have about 49% IR character and very large oscil- 
lator strength of 5'a=3.89. However, this large Raman- 
IR mixing might come from that this modes accidentally 
degenerate with the 407 and 433 cm~^ IR modes in the 
numerical calculations, and might not be the case in the 
real system. Unfortunately, we can not find suitable ex- 
perimental IR spectra to compare with for the Bi modes. 
We hope future experiment can clarify this problem. 



2. Raman spectra 

The temperature-dependent Raman spectra of 
HoMn205 and TbMn2 05 were measured by Mihailova 
et al.^^ We compare the calculated and measured Raman 
spectra of TbMuaOs in Table IVH El |VlT] and IVlIIl for 
the Big, Ag, B^g and B2g modes respectively. The A^ 
silent modes are also given in Table IVIIIl for they also 
Raman-active in structure L. Some Raman modes were 
missing in the experiments probably because the missing 
modes are of very low intensity or out of the spectral 
range {lu < 100c?7i^^) of the experimental setupi^ Once 
again, the calculated phonon frequencies are in very 
good agreement with the experimental values. 



TABLE VIII: Calculated phonon frequencies uj of A2 modes, 
classified according to their major symmetry characters and 
compared with the experimental data. The Raman-active 
phonons are extracted from Ref. |24| . 



UJ 

GGA 


(Silent) 
Raman% 


B 

UJ 

GGA 


2g (Raman) 

Raman% 
Exp. 


110.2 
131.4 


0.04 
0.09 


95.9 
210.6 


214 


99.9 
99.4 


226.6 
293.9 


6.18 
40.7 


231.1 
274.2 


232 


99.9 
99.4 


399.4 


12.6 


288.9 


301 


60.3 


425.0 
497.7 


3.00 
5.05 


337.5 
441.1 




99.3 
89.1 


529.4 
611.5 


1.96 
5.65 


448.7 
469.6 


455 
470 


96.9 
92.9 






478.6 
562.6 


505 


99.3 
99.1 



So far, no IR-active modes has been observed in the 
Raman spectra of TbMn205 . However, there is an ev- 
idence that one of the IR modes become also Raman- 
active in HoMn205j2i Three low frequency Ag Raman 
modes were observed at 10 K in HoMn205r^ and they 
are the 217 cm^^, 219 cm~^ and 226 cm~^ modes, 
whereas the 219 cm^^ mode disappears at room tem- 
perature. Since the phonons of frequencies at this range 
are mainly due to the vibrations of Mn and O, and less 
related to the R atoms, one could therefore expect there 
are also three Raman modes in this frequencies range for 
TbMn205. However, only two Ag modes (215 cm~^ and 
221 cm— 1 modes) were observed (nevertheless at T— 300 
K) in this frequency range in TbMn205, consistent with 
our first-principles calculations for the Ag modes. We 
therefore propose that the 219 cm~^ mode in HoMn205 
is actually an IR-active phonon that become also Raman- 
active at low temperature. 



VI. DISCUSSIONS ON THE MICROSCOPIC 
ORIGIN OF THE ELECTRIC POLARIZATIONS 
AND GIANT MAGNETOELECTRIC COUPLING 

We have now firmly established that the ground state 
structure of TbMn205 {L or R) is indeed of symmetry 
Pb2im, allowing spontaneous polarization along the b- 
axis. We have calculated^- the spontaneous polarization 
of structure L (R) via Berry phase techniquci^^ The in- 
trinsic polarization in this material is calculated to be 
1187 nC-cm~^." This value is an order of magnitude 
smaller than that of the traditional ferroelectrics-l (e.g, 
for BaTiOs ~ 25/iC-cm~^), yet it is about 30 times larger 
than the currently measured valued {^ 40 nC-cm^^ ) for 
this compound. The reasons for the discrepancy between 
the calculations and experiments might because of the 
approximations used in calculations. For example, we ap- 
proximate the magnetic propagation vector kz=0.32 by 
zero. We also ignore the spin-orbit coupling and assum- 
ing coUinear spins. Without these approximations, the 
polarization might be reduced. On the other hand, the 
experiment'^ measured polycrystalline samples, in which 
the grains polarize along different directions cancel each 
other, therefore might significantly underestimate the in- 
trinsic polarization. We believe a high quality single crys- 
tal sample should greatly enhance the measured electric 
polarization. 

To further elucidate the origin of the polarization, we 
calculate the spontaneous polarization for the high sym- 
metry structure H and get P= 228 nC-cm^^. It might be 
surprising to see a crystal that possesses inversion sym- 
metry would develop spontaneous polarization at the first 
sight. However, to discuss the symmetry of a crystal with 
a magnetic ordering, the pure point groups may not be 
enough. One has to resort to the magnetic groups to take 
account of the spin components. Based on the symme- 
try considerations, there are three possibilities (but not 
exclusive) to develop spontaneous polarization: 

(i) The crystal structure itself does not have inversion 
symmetry. This is the most familiar case to us, as it is 
seen in the traditional ferroelectrics. 

(ii) The crystal structure itself has the R~^ symmetry, 
but the magnetic structure does not have the R~^ sym- 
metry and the combination of the R^^ and time- reversal 
(r~^) symmetry, (RT)^^. If the magnetic structure of 
a system does not have the R~^ symmetry, but still has 
the {RT)~^ symmetry, each spin channel may have finite 
electric polarization. However, the polarizations from 
spin up and spin down channels cancel each other, lead- 
ing to zero total electric polarization. 

(iii) Even both the crystal and magnetic structure have 
the inversion symmetry, but they do not have the com- 
mon inversion center, as proposed by Betouras et al.^ 



In the later two cases, the "electronic" inversion sym- 
metry is broken because of the magnetic ordering. The 
electronic symmetry breaking will further couple to the 
lattice and lead to lattice distortion. It was shown^^ that 
in TbMn2 05, the SC g does not have the R^^ and the 
{RT)~^ symmetry (case ii), because the special spin con- 
figurations of the Mn'^+-Mn^+-Mn^+ chains along the b 
axis. Therefore, when holding atoms fixed at the cen- 
trosymmetric structure, turning on the magnetic order 
does two things: it generates a purely electronic polar- 
ization of 228 nC-cm~^, and it also applies forces to the 
atoms. These forces give rise to atomic displacements 
that yield an additional 959 nC-cm~^. 

The above discussion is based solely on the symmetry 
argument. We now discuss the microscopic mechanism 
of coupling between the magnetic order and lattice. The 
microscopic origin of the ME coupling may come from^ 

(i) Symmetric superexchange interactions that of the 
scalar field type: Sn • Sn+i. 

(ii) Antisymmetric superexchange interactions that of the 
wector field type: Sn x Sn+i. 

The antisymmetric superexchange interaction (also 
known as Dzyaloshinskii-Moriya interaction) originate 
from the spin-orbit interaction. For the antisymmetric 
superexchange mechanism, nonocolhnearity is essential 
for the ME coupling. Experimentally, it was shown that 
in TbMn205 the largest electric polarization is associ- 
ated with a commensurate magnetic phase that is almost 
coUinear ^ii^ therefore, in this phase, the spin-orbit con- 
tribution to the ME coupling is small. An alternative 
mechanism must be considered to account for the ME ef- 
fects. In the present calculations, we ignore the spin-orbit 
coupling and assume coUinear spin (S^ || Sj), therefore 
the symmetry breaking and ME coupling is merely due 
to the symmetric superexchange Sn • S„+i. 

We now derive an effective model to describe the mi- 
croscopic mechanism of the ME coupling in TbMn205, 
based on a Heisenberg model assuming NN interac- 
tions. The total energy of the system around the high- 
symmetry structure H can be written as, 



Ei{u,}) = Eo + lj2 



mxuj\u\ 



A i'j 

(8) 
where ux is the A-th zone-center phonon, and S^ is the 
magnetic moment of the i-th atom. Here, we consider 
only the magnetic moments of the Mn'^+ and Mn^+ ions. 
Eq is the total energy of structure H without spin-spin 
interactions. Since the magnetic-ordering-induced lattice 
distortion from structure H is small, we could expand the 
exchange interaction Jij{{ux\) up to the second order of 
Ux, i.e.. 
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E{{ux}) = {Eo-J2Jv(0)S,-Sj 



EES«s..s, 



r 



:^E™> 



2 2 



EE 

!j Xp 



du\du 



-uxUpSi-Sj). (9) 



The above three terms play different roles in the multi- 
ferroics. 

(i) The zeroth order term in u\, 



Eq — 2_^ JijiO)^i ■ Sj , 



(10) 



determines the ground state spin configurations which is 
further discussed in Appendix [K\ 



(ii) The second term in Eq. p4p that is linear in u\, 

(11) 



EE?^-^s..s,, 



du). 



provides the driven forces for the lattice distortion from 
high symmetry structure H. The force on the A-th 
phonon modes is 



dE 



EuJij 



'7^i • Oj 



(12) 



which is nonzero, provided {8^} does not have the R~^ 
and {RT)~^ symmetry. 

In Ref. [l6|, we illustrated the lattice distortion from the 
high symmetry structure using a spin chain model. For 
SC g, inside the Mn^+-Mn'*+-Mn^+ segments along the 
6-axis, Mn^"*" always have the same spin with the upper 
Mn'^+ but opposite spin with the lower Mn^+, whereas 
for g', the opposite is true. Mn''+ could therefore move 
closer to the Mn'^+ with same spin to lower the exchange 
energy. According to Eq. p2|) . wehaveFx(g) — —F\{g'). 
The high symmetry structure H then spontaneous relax 
to structure L ov R according to its SC. The lattice dis- 
tortion from high symmetry structure can be estimated 
as, 



A ^J 1-2 



(13) 



where uj is the weight averaged phonon frequency. We 
could therefore expect that a material with larger dJ/du 
and softer phonon frequencies would have larger ME ef- 
fects, provided the symmetry requirements are satisfied, 
(iii) The quadratic term in u\, 



1 



"^^ mxujlul 



2^^ 

X 



EE 

ij Xp 



-^J^uxup^, ■ S, (14) 



describes the coupling between phonons and magnons, 
which renormalizes the phonon frequencies. For exam- 
ple, in ZnCr204,— i^ a large splitting between IR ac- 
tive phonons along the x and z axis was observed due 



to the phonon-magnon coupling. The phonon-magnon 
coupling might be responsible to the "step" in dielectric 
constant at 24K in TbMn205 . Detailed investigation 
of the phonon-magnon coupling in TbMn2 05 is out of 
the scope of the present work, and we leave it for future 
investigations. 



VII. SUMMARY 

We have investigated the ground state structural and 
electronic properties of multiferroic TbMn205 using den- 
sity functional theory within GGA approximations to 
elucidate the microscopic origin of the ferroelectrcity and 
giant magnetoelectric coupling. We use the coUinear spin 
approximation and ignore spin-orbit coupling. The calcu- 
lated lattice and electronic structures agree well with the 
known experiments. We confirm that the ground state 
structure of TbMn205 is of space group Pb2im, allow- 
ing polarizations. The spontaneous electric polarization 
is calculated to be 1187 nC-cm"^ along the b axis. The 
calculated zone-center optical phonons frequencies and 
the oscillator strengths of IR phonons agree very well 
with the experimental values. Unlike the traditional fer- 
roelectrics, no soft phonons is found in TbMn2 05. We de- 
rive an effective Hamiltonian to explain the magnetically- 
induced ferroelectricity, in which the spin-lattice cou- 
pling is of symmetric superexchange interaction type. 
Our results strongly suggest that the ferroelectricity in 
TbMn205 is driven by the magnetic ordering that breaks 
the the inversion symmetry without invoking the spin- 
orbit coupling. 
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APPENDIX A: GROUND STATE SPIN 
CONFIGURATIONS 

In this appendix, we discuss in more details about 
the magnetic structures of TbMn205 and the (su- 
per)exchange interactions Jsi^ To determine the Js, we 
calculate the total energies of eight SCs shown in Fig. 
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TABLE IX: The DFT-GGA calculated total energies of spin 
configurations (shown in Fig. [5} under the high symmetry 
structure H, where E3 = — 2 J3 S3 • S4,, E4 = — 2 J4 S3 ■ S4, 
E5 — — 2J3 S3 •S3. The exchange integrals J3, J4 and J5 are 
fitted using the total energies of four AFM structures. AE 
are the energy differences between the fitted energies from the 
Heisenberg model and the calculated ones. Unit of energy is 
in meV, and the unit of exchange integrals are in laeV/ 11% . 



sc 


status 






Energy 


AE 


a 


FM 


Eo + I6E3 + IGE4 + 4E5 





-179.84 


b 


FM 


Eo- 


16Ei - 16^4 + 4E5 


-1355.66 


-158.22 


c 


FM 




Eo - 4E5 


-1004.93 


-8.41 


d 


FM 




Eo 


-910.32 


-19.78 


e 


AFM 




Eo + 8E:i + 8E4 


-596.59 





f 


AFM 




E(] — SE's — 8-E'4 


-1263.61 





a 


AFM 




Eo - 16^4 - 4^5 


-1624.95 





h 


AFM 




Eo + 16^4 - 4^5 


-407.72 









J3-- 


=-0.45 ./4 = -4.92 


Js = -1.85 





[51 under the fixed high symmetry structure H. The to- 
tal energies of these SCs are listed in Table IIXI among 
which SC g has the lowest total energy. We then use a 
Heisenberg model with NN interaction, 



E ~ Eo 



/ , "Jij '^i 



(Al) 



to fit Js and to check whether the g is indeed the the 
ground state SC. Here, Si is the spin of the i-th Mn 
ion and Jij denotes the exchange integral between two 
NN i-th and j'-th atoms. J3, J4 and J5 are fitted using 
four AFM SCs, and the results are listed in Table [IXl 
Since Ji and J2 is not relevant here, they are not fitted. 
We calculate the total energies of the four FM SCs from 
Eg. ljAip . These results are in reasonable good agreement 



with the first-principles calculations, suggesting that the 
Heisenberg model is valid even in the FM phase when 
applied to TbMnsOs. From Table [IXl we find all the 
three exchange integrals are negative, i.e., of AFM type, 
and I J3I ^ IJ4I, I J5I. Therefore the ground state SC are 
mainly determined by J4 and J5. It is easy to see that 
SC g (and g') indeed has the lowest energy in this model. 

In both SCs g and g' , Mn'^^ form an AFM square 
lattice in the ah plane, whereas Mn'^^ couples to Mn'^^ 
either antiferromagnetically via J4 along the a axis or 
with alternating sign via J3 along the h axis. The zigzag 
chains along the h axis can be simplified as Mn'^+-Mn''+- 
Mn'^^ segments linked via J5 superexchange interactions. 
Inside the segments Mn'^"'" and Mn'*+ interact through su- 
perexchange J3, whereas the chains couple to each other 
through J4 . The AFM chains along the a axis are of two 
types, labeled as I and II respectively in Fig. [3] . Adja- 
cent chains I and II along the h axis couple to each other 
via J3, in which half of them are frustrated. 

The SCs g and g' have several energetically degenerate 
SCs as shown in Fig. [31 There are two relative positions 
between chain I and II, / and r. "/" ("r") means that the 
Mn'^^ ions in chain I have the same spins to the Mn'*+ 
ions on their left [right) side in chain II in each unit 
cell. Figure [31 shows 4 different combinations of chains 
I and II. Figure [3ja) and (b) are the ground SC g and 
g' respectively, whereas Fig[3fc) and (d) are in fact two 
domain wall structures dw and dw' . The two domain 
walls are energetically degenerate to g and g' in the NN 
Heisenberg model, because the J3 exchange interactions 
between chain I and II cancel out no matter what their 
relative positions are. However, after turning on the spin- 
lattice coupling, the crystal structure H will relax to L or 
R, and the energy degeneracy between the domain wall 
structures {dw and dw') to SCs g and g' will be lifted. 
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